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TABLE OF SYMBOLS 



Below are two tables of symbols. The first table lists 
all the symbols which occur in the text and, where applicable, 
their FORTRAN counterparts used in the EULER-1 code. The 
second table lists the remaining symbols used in the EULER-1 
code which were not listed in the first table. Both tables 
are in alphabetical order. 



SYMBOLS USED IN THE TEXT 



Text EULER-1 

A A 

P PRESS 



Q QQ 

q Q 

q r 

R RR 



S 



S 



s, n,m 

T TEMP 

t T 

u U 



V VS 

s 

W w 



Definition 
Speed of sound 
Static pressure 

Modified Riemann variable (Q = q +AS) 
Absolute velocity magnitude 
Reversible heat transferred 
Modified Riemann variable (R = q -AS) 

Gas constant 

A modified form of entropy 

Unit vectors in the s, n, and m directions 

Static Temperature 

Time 

Velocity magnitude relative to a steady 
shock wave 

Shock wave velocity 

Mach number relative to a steady shock 
wave 
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w 

Z 

X 

A 

6 



e,<j> 



p 



Vector of the principle variables, Q,R,S 

Z (K) Vector of the right hand sides of the 

governing equations 

LMD(K) The characteristic trajectories in the 

space-time plane (q+A, q-A, q) 

Small spatial change 

Small change with respect to time 

Flow angles with respect to reference 
coordinate planes 

DENS Density 



Y G 



Ratio of specific heats 



Symbol 

*A 

ABAR(K) 

*B 

AR 

COUNT 
D ARRAY 

DEL* *H 

DELX (K) 

DEI,**L 



ADDITIONAL SYMBOLS USED IN EULER- 1 

Definition 



Suffix which denotes the left side of a 
discontinuity (* may be any variable) 

The average value of A between the limits of 
integration of Z (K) 

Suffix which denotes the right side of a 
discontinuity (* may be any variable) 

The ratio of the sound speeds across a shock, 
(high/low) 

A counter for the number of time steps 

The array of density values to be plotted by 
the graphics routines 

The change in the variable ** in going from node 
I to node 1+1 [ **(I+L) - **(I) ] 

t - h 

The distance from the point where the K 
characteristic crosses the known time level to 
the I-h node, measured positive to the right 

The change in the variable ** in going from node 
1-1 to node I [ **(T) - **(I-1) ] 
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DLCD 


The exact value of the density to the left of a 
right propagating contact discontinuity in the 
Riemann problem 


DLSH 


The exact value of the density to the left of a 
right propagating shock in the Riemann problem 


DLTA* * 


A prefix which indicates the spatial change in 
** over one time step 


DQ 


The change in the velocity magnitude across a 
shock, (high-low) 


DR 


The ratio of the density across a shock, (low/high) 


EE 


A value which indicates the precision to which the 
characteristics are to be calculated 


E(K) 


The computed error in the calculation of the 
characteristics 


GRAPHS 


Parameter which controls the type of output pro- 
duced by EULER- 1, 0 = Tabular, 1 = Plot of Q, S, 

pressure and density, 2 = Comparison of density 
with exact solution 


G1 


1/ (G-l) 


G2 


2/ (G-l ) 


H 


The non-dimensional value of one spatial interval 
1/ (N-l ) 


I 


Subscript which denotes the spatial node 


* * I NT ( K ) 


The interpolated value of ** at the point on the 
known time level where the K^h characteristic 
crosses 


INTEG (K) 


The result of integrating Z (K) 


12 (L) 


The node to the right of discontinuity L 


J 


The time level 


JSTOP 


The number of time levels desired to be computed 


K 


A subscript which indicates the characteristic 
being dea] t with. 1 = Q+A 2 = Q-A 3 = Q 
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L 


A subscript which denotes the type of discontinuity 
1 = schock 2 = contact discontinuity 3 = head 
of rarefraction wave 4 = tail or rarefraction 
wave 


M 


Parameter which controls the types of discontinui- 
ties which EULER-1 tracks. 1 = shocks only. 

2 = shocks and contact discontinuities 


N 


The number of nodes in the computational grid 


NEW** 


A temporary storage location for variables ** 
updated in time 


PARRAY 


The array of densities plotted by the graphics 
routines 


PR 


The ratio of pressures across a shock, (high/low) 


PRI 


The initial pressure ratio across the diaphragm 
(high/low) 


*PRIM (K) 


Suffix which indicates the spatial derivative of 
* at the present time level, where * is either A 
or Q 


QARRAY 


The array of velocities to be plotted by the 
graphics routines 


QLI 


The initial value of the velocity to the left of 
the diaphragm 


QRI 


The initial value of the velocity to the right of 
the diaphragm 


QQJO 


The measured change in the Riemann variable QQ 
across a shock (high-low) 


QQJE 


The change in the Riemann variable QQ across a 
shock calculated analytically (high-low) 


SIGMA (L , J) 


Discontinuity locations. J = 1 indicates the 
known time level and J = 2 indicates the unknown 
time level 


* *STEP 


The change in the primary variables with respect 
to time during one time step (** = QQ, RR, S) 
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SARRAY 


The array of entropy values to be plotted by the 
graphics routines 


SKIP 


The number of time steps between outputs 


TRI 


The initial temperature ratio across the diaphragm 
(high/low) 


VCDE 


The exact velocity of the contact discontinuity 


VHEAD 


The exact velocity of the head of the rarefraction 
wave 


VTAIL 


The exact velocity of the tail of the rarefraction 
wave 


VS 


The shock wave velocity computed by EULER-1 


VSE 


The exact velocity of the shock wave 


X 


The non-dimensional spatial position 


XARRAY 


The array of node locations used by the graphics 
routines 


XINIT 


The initial location of the discontinuities. Used 
to calculate the exact solution 


X2 (L) 


The position of the node to the right of a 
discontinuity 



11 



ACKNOWLEDGMENT 



Thanks to my wife, Kathy, for her support and encourage- 
ment, to Augie Verhoff, for his personal assistance, to Atul 
Mathur, for his availability and assistance, and to Ray 
Shreeve, for making this an enjoyable and worthwhile learn- 
ing experience. 



12 



INTRODUCTION 



I , 

A program is underway at the Naval Postgraduate School's 
(NPS) Turbopropulsion Laboratory to evaluate the wave rotor 
concept. The wave rotor, operating as a component in a gas 
turbine engine, uses unsteady wave propagation in tube-like 
passages to compress incoming air before it goes to the com- 
bustion chamber. The combustion chamber output is routed back 
to the rotor to create the unsteady waves. Thus the rather 
simple rotor, with "partial admission" inlet and outlet ports, 
acts both as a compressor and as a turbine. The alternation 
of hot and cool gases through the same hardware aids cooling 
and allows higher operating (cycle) temperatures to be used. 

The interaction of the hot and cool gases within the passages 
of the wave rotor through the wave patterns that are created 
pose a difficult problem. The work underway aims to develop 
and validate preliminary design and performance analysis tools. 

One of the tools needed is a computer program which can 
be used to construct wave interactions in a design process 
and then accurately predict the performance of the design. 

Until efficient and accurate methods of solving the full 
Navi er-Stokes equations for unsteady turbulent flows are 
developed, the design program will be based initially on the 
solution of the unsteady Euler equations and a method devised 
to represent losses. Once the program is developed, it must 
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be verified against experiment. This is the overall goal of 
the present program in which a wave rotor apparatus has been 
assembled and methods of measuring the unsteady pressures and 
temperatures are being developed concurrently with the compu- 
tational effort. 

Three different approaches to the solution of the unsteady 
Euler equations were examined in the overall program. First, 
Eidelman developed a tv/o-dimens ional code based on the 
Godunov method of solution [Ref. 1] and applied the code 
to examine unsteady wave propagation in ducts [Ref. 2] and 
the process of port opening to wave rotor passages [Ref. 3] . A 
summary of Eidelman' s work is given in Reference 4. While 
the method is conservative and does not require the introduc- 
tion of artificial viscosity, the extension from one to two 
dimensions is not rigorous when shock waves are present, and 
computational times with the present code are quite long. 

The extension to include viscous effects would require a separate 
treatment of the boundary layer. 

Second, Mathur developed a one-dimensional code based on 
the Random Choice Method (RCM) of solution [Ref. 5]. Somewhat 
similar to the Godunov method, in that the solution is based 
on solving the Riemann problem within each grid cell at each 
time step, the RCM approach results in very sharp discontinui- 
ties which can be tracked easily. This is particularly useful 
in constructing wave rotor cycles, in which the position of 
the gas-gas interface is equally as important as the position 
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of the compression and expansion waves. The code is there- 
fore valuable in the preliminary design process, to examine 
suitable port arrangements, and the gas properties at the 
ports, for a given task. Unfortunately, an extension to two- 
dimensional and/or viscous flow can not be made rigorously. 

The third approach was followed in the present work. The 
QAZ1D method for compressible inviscid flow computations 
developed by Verhoff and O'Neil [Ref. 6] was implemented to 
generate a one-dimensional unsteady Euler code with the goal 
of evaluating the suitability of the approach for wave rotor 
applications. Some advantages of the QAZ1D method were recog- 
nized as being the following: 

1. The method is based on the use of characteristics. 

Such methods can model wave propagation accurately. 

2. The use of a natural streamline coordinate system 
eases the difficult task of computing with two and 
three dimensional grids. 

3. The equations are written in a form which allows a 
straightforward extension to viscous flows. 

4. Codes for computing .internal, steady, two-dimensional 
flows in the presence of shocks, and simple internal 
viscous flows, have been generated quickly without 
significant development problems [Ref. 6]. 

In the work reported herein, a one-dimensional FORTRAN 
code based on the QAZ1D method was developed using the NPS 
IBM370-3033 computer and subsequently exercised on the shock 
tube test problem. In reporting the work, first, in Section 
II and Appendix A, a complete account is given of the deri- 
vation and non-dimensionalization of the governing equations. 
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In Section III, a FORTRAN code, EULER-1, is described. Its 
operation on the NPS computer, and the listing of the code, 
are given in Appendix B and Appendix C, respectively. 

Results of applying the code to the shock tube test problem 
are given in Section IV. Difficulties encountered in the 
implementation of the method and additional comments are given 
in Section V, and conclusions are given in Section VI. 



16 



II. THE QAZlD METHOD 



A. OVERVIEW 

The QAZlD method uses Riemann-like variables with a modi- 
fied entropy term in their definitions, to express the Euler 
equations in a natural streamline coordinate system. The 
resulting partial differential equations (PDE) , when cast 
along the characteristic trajectories in the space-time 
domain, reduce to a system of ordinary differential equations 
(ODE) which may be solved explicitly. The advantage in using 
the modified Riemann variables is that they are less affected 
by discontinuities in the flow than are the standard Riemann 
variables. Since the equations are not valid across discon- 
tinuities which cause irreversible losses, such discontinui- 
ties must be located and treated with special logic in the 
numerical formulation. 

In the following presentation, the reader's familiarity 
with the method of characteristics is assumed. 

B. DEVELOPMENT OF THE EQUATIONS 

This section presents the governing equations and outlines 
their development. A rigorous derivation of the equations 
is presented in Appendix A. 

1 . The Coordinate System 

The natural streamline coordinate system (s,n,m), is 
shown in Fig. A-l relative to a fixed rectangular cartesian 
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system (x,y,z). The (s,n,m) system is a right-hand orthogonal 
system which undergoes curvilinear translation as it moves 
with a fluid particle along a streamline. The coordinate 
system is described in detail in Appendix A. 

2 . Variables 

The modified Riemann variables, or "extended Riemann 
variables" [Ref. 6:p. 1], are defined as 



Q = q + AS 
R = q - AS 



( 1 ) 



where q is the velocity magnitude, A is the speed of sound 
and S is the modified entropy defined in terms of pressure 
and density as 



S 




Y (Y-l) 



[ 2y - £n(P/p y )] 



( 2 ) 



The modified entropy relation is the result of defining the 
entropy change, dS , to be given by 



dS = 




( 3 ) 



where dQ^ is the heat required to be added in a reversible 
process between the same end states, T is the temperature and 
Y is the ratio of specific heats. 
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3. 



Conservation of Mass 



By applying the continuity equation to a differential 
stream tube of variable cross sectional area, in natural, 
coordinates, and ignoring all third order and higher deriva- 
tives, one obtains 



3p 

3t 



+ q 



3s 



3q 

3s 



r 3 9 Q 

p q [ + cos 0 



3 d 

3m 



] = 0 



(4) 



where q is the velocity, p is the density, s is the stream- 
wise spatial dimension, and 0 and <j> are the flow angles as 
defined in Fig. A-l of Appendix A. 

4 . Conservation of Momentum 

Applying the vector form of the momentum conservation 
law in natural coordinates, the equation of motion for inviscid 
flow becomes 



? r 9q . 3q 3P, 

V p at + p + a ! 1 



? , 3 0 2 30, 3P-, 

+ L n [p q 3t + p q 37 + 3^ ] 



(5! 



? r Q ii . 2 ^ 3j> , 3P, 

+ I m tp q cos BjJ+pq cos 6 3t + 3i^ ] 



= 0 



5 . The Conservation of Energy 

In the absence of friction and heat conduction, and 
outside of discontinuities (through which irreversible 
changes consistent with mass, momentum and energy conservation 
are permitted) , energy conservation is equivalent to the 
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statement that the entropy of a fluid particle does not change. 
In the natural coordinate system, therefore 



dS _ as as 
dt at q as 



0 



( 6 ) 



Equation (6) will not be valid across shock waves or contact 
surfaces between gases having different states. 

6 . Transformation to a Useful Form 

Equations (4) and (5) are first expressed in terms 
of the primary variables q, A, S, and the flow angles. Using 
the definition of sound speed in a perfect gas 

A 2 = YP/p (7) 



Eq. (4) becomes 



3A 

at 



+ 



aA 
q al 



Y-l 



A la 

A as 



Alzlrll 

A 2 l 3n 



+ cos 9 



Mi 

3m J 



- 0 



( 8 ) 



With Eq. (3), the equation of state for a perfect gas,, and 
the first law of thermodynamics, Eq. (5) gives 



3q 

at 



+ 



q 



la , _2A_ aA 

as y - i as 



2 as 

1 a¥ 



o 



ae 

at 



+ 



q 



30 

as 



a<ft 

at 



+ 




A 2 3 £n P 

Yq an 

A 2 a f.n P 
Y q cos 6 am 



( 9 ) 
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Equations (6), (8) and (9) constitute the system of P.D.E.'s 

which describe the unsteady isentropic flow of a perfect gas 
in natural coordinates. Since third order terms were 
neglected in the conservation of mass, the system is only 
accurate to second order. 

The system of P.D.E.'s is now transformed into a system 
of O.D.E.'s along the characteristic trajectories. In general, 
if w is a function of (s,t), then 



dw 



3w , , 3w 

35 ds + at dt 



or 



dw _ Sw.ds, _3w 

dt 3s dr 3t 



( 10 ) 



where (ds/dt) describes the "characteristic" direction in the 
(s,t) plane along which Eq. (10) gives the rate of change of 
the parameter w. After the appropriate algebra and the intro- 
duction of Eq. (1), the final system of equations becomes 



f + (< 3 +A) !§ = 



1 - a (s ?_) ^ ) i 

2 A lb y-1 ' 3s ' q y-1 A } J 



^ q A s iff + c °s e f|] 



3R , . . 3R 

3t + (< 3" A) 35 



= + 



y-1 



2 . r 3 



A(s -Fi )1 ^ (q+ 5h: A!1 



+ 1 i iqA s l ls * cos 6 l^ ] 






ID 
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3S 
3 1 



+ 




0 



i! 

3 1 



+ 



q 



ii 

3s 



A^_ 3 Jin P 
yq 3n 



ii + 

3t 



- 



A 2 3 In P 
y q cos6 3m 



J 



The characteristic directions in the space-time plane are 
clearly q+A , q-A , and q. 



C. 



SOLUTION METHOD 

Equation (11) may be expressed as 



3w 

3t 




Z 



( 12 ) 



or along the A directions as 



dw 

dt 



Z 



(13) 



where the vectors are defined as 



r i 

Q 




— — 

q+A 




R 




q-A 




< s 


II 


q 


Z = right hand 


0 




a 


side of 






i 


Eq. (11) 


■ 0 , 




q 
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A trajectory in the space-time domain is illustrated in Fig. 1, 




► 6 1 

[6w = w fi -w ] 
i 



J 



l Aw = w — w ] 
s i A 



Figure 1. Solution Procedure in the Space-Time Domain 

which shows an infinitesimal interval of space between two 
"nodes" of the computational mesh. 

For the change in w from A to B .along the trajectory 
with slope X, 



t + 5t _ 

/ Z dt 




w 



A 



[w^ - w ] + [w - ] 

L B s . s . A 

l l 



- 5w + Aw 
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where 6w denotes the change due to time at a fixed location 
and Aw denotes the change due to displacement at a fixed 
time, for the characteristic trajectory. The essence of the 
solution procedure is to calculate the change in the varia- 
bles at each spatial node during the differential time 
interval, so that the step to the next time level can be 
made. In other words, 6w must be calculated at each node 
using 



_ t+5t _ 

6w = -Aw + / Z dt (14) 

t 

Since all the information at time level t is known, the 
position of A can be determined by an iterative procedure 
based on X and Aw can be calculated by interpolation. The 
line integral can be evaluated by transforming the integral 
to a purely spatial integral using X = ds/dt. Thus 

t+6t _ B - 

Z dt = It ds = 
t A s^+As 

and the integral can be evaluated by any one of several 
numerical methods. Thus, the variables at each node can be 
updated in time using Eq. (14) and the process repeated. 
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D. DISCONTINUITIES IN THE FLOW 

There are several types of discontinuities that must be 
considered. They are gradient discontinuities, contact 
discontinuities and shock waves. 

Gradient discontinuities are characteristic of the head 
and tail of rarefraction waves, the collision of two shocks and 
the interaction of a shock and a contact discontinuity. 

Across such a gradient, the derivatives of velocity, pressure 
and sound speed are discontinuous. 

A contact discontinuity is caused by the interaction of 
two shocks of opposite family and when a shock overtakes a 
shock of the same family. Entropy and sound speed are discon- 
tinuous at a contact discontinuity. 

Across shock waves, velocity, pressure, density and 
entropy are discontinuous [Ref. 7], 

Since Eq. (11) is not valid across discontinuities, addi- 
tional logic is required in the numerical procedure to model 
flows which contain discontinuities. The method presented in 
Reference 6 for making a correction in the case of shock waves 
will give good accuracy in problems where the solution is 
converging to a steady state condition but will not give 
accurate results during the transient portion of such problems 
or for problems with only unsteady solutions. In the case of 
applications to the wave rotor, it would certainly be necessary 
to know the locations of the contact surfaces between the 
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hot and cool gases as well as the locations of the shock and 
expansion waves. 

The methods used in the present work to correct for dis- 
continuities are those of Moretti [Ref. 7], The present 
discussion will be limited to the treatment of shock waves. 
The method makes use of the analytical relationship between 
the change in the Riemann variables across a shock and the 
incoming Mach number relative to the shock wave (W) which is 
illustrated in Figs. 2 and 3. The relation is used to 
determine the shock speed and to transform the problem to 
a steady case which can be handled rising normal shock rela- 
tions. For the situation depicted in Fig. 2 of a shock 
propagating to the right with velocity V g into air with 
velocity q fi , and with the high pressure side to the left, 
where A denotes the left side of the shock and B the right, 
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If the pressure and density are non-dimensionalized by the 
values on the low pressure side of the shock, the change in 
the extended Riemann variable Q is given by 
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where 



4 = 



1/2 



This exact relationship can be approximated by the polynomial 



Q A _Q B 2 

, = -2.7574 + 3.1573W - 0.2863W 

A B 



(i7: 



as illustrated in Fig. 4 over a shock strength range of 1.0 
to 4 . 0 . It should be noted that if the high pressure side 
were to the right, as in Fig. 3, Eq. (15) would become 



W = 



U A 

*1 



q A -v 

A s 



and the left side of Eq. (16) would become 



R a - R n 
A B 



The procedure then, is to measure the change in Q across an 
interval where a shock is known to exist. This value, call 
it AQ^, is used in Eq. (17) to get an approximation for the 
corresponding value of W and then the exact value of AQ, call 
it AQ^, is calculated using Eq . (16). If the exact value of 

AQ is not equal to the value measured across the shock, a 
new value, AQ, is calculated according to 
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AQ. + (A Q -AQ) 
i m E 



and entered into Eq. (17) to obtain another value of W. When 
the exact value of AQ calculated from Eq. (16) equals the 
measured value of AQ, W is known to be correct and the normal 
shock relations can be used to calculate the values at node A. 
Also, V g can be calculated from Eq . (15) and used to track 

the shock during the next time interval. 
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III. 



DESCRIPTION OF FORTRAN PROGRAM " EULER- 1" 



A. MACHINE AND LANGUAGE 

The EULER-1 program is written in VS FORTRAN and runs on 
an IBM 3033, System 370 computer, however the code is simple 
and small enough to enter and run on a mini or micro computer 
in a Basic language if one is willing to accept significantly 
longer run times. Table 1 at the end of this section is a 
summary of the editable parameters and their effect on the 
program. 

B. CONVENTIONS AND BASIC STRUCTURE 

EULER-1 is a first-order one-dimensional code using an 
evenly spaced numerical grid. All values are double precision 
except those used in the graphics routines which are rounded 
to single precision. 

Each subroutine has its own variables space. In other 
words, variables are not shared in common throughout the 
program but must be passed to the subroutine being called 
by the calling routine. Tn all cases, however, the variables 
have the same name in both the called and the calling routine 
so there should be no confusion. This was done so that 
arrays could be dimensioned at execution time. 

The program, depicted in Fig. 5, is structured around a 
main routine which serves as a user input area, sets up the 
problem and calls five subroutines to solve the problem and 
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a) Main Program 



Figure 5. EULER-1 Flowchart 
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b) "SWEEP” Subroutine 
Figure 5. (Continued) 
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output the solution. The five subroutines are. "TIME," "TRAK, " 
"SWEEP," "JUMP" and then one of several output routines depend- 
ing on user desires. Output options include two types of 
graphical displays, "PLOT" and "EXACT," and one tabular 
listing routine, "LIST." When "PLOT" is selected, "BORDER" 
is automatically called to set up the plotting area. There are 
two other subroutines, "RSHOCK" and "LSHOCK," which are called 
by "SWEEP" as needed. 

C. SUBROUTINE DESCRIPTIONS 

In general, each subroutine begins with a heading followed 
by variable definitions where appropriate, variable declara- 
tions and array dimensioning. Most variables are defined in 
the "MAIN" routine and only those variables which were not 
are defined in the subroutines where they are used. 

1 . The "MAIN" Routine 

This routine forms the main structure of the program 
and includes a heading, an extensive list of variable defini- 
tions, a. user input area, the necessary statements to make 
the initial value assignments to variables, and the call 
statements for the various subroutines. 

The input area is those lines of the program (145-175) 
where the user edits the program to establish the initial 
conditions, mesh size, termination criteria and output options. 

The initial conditions which may be modified are the 
pressure, temperature and density ratios across the diaphragm, 
the initial velocity of the fluid on each side of the 
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diaphragm, and the value of gamma (which must be the same for 
both sides) . 

All velocities are non-dimensionalized by the initial 
sound speed on the low pressure side of the diaphragm and 
pressures and densities are non-dimensionalized by their 
initial values on the low pressure side of the diaphragm. 

The mesh size may be set at any odd number and the 
arrays in the user input area must be dimensioned as such. 

The criteria for program termination is the number of 
time steps computed, which the user selects. 

The output options are determined by the value of the 
variable GRAPHS. A value of zero results in a call to "LIST" 
which writes to file 9 on the. .user's permanent disk, a tabular 
listing of the variable arrays and discontinuity locations. 

A value of 1 results in a call to "BORDER" and "PLOT" which 
create a plot of the pressure, density, velocity and entropy 
distributions as in Fig. 7. A value of 2 results in a call 
to "EXACT" which creates a plot of the density distribution 
compared to the exact solution as in Fig. 8. The exact values 
to be plotted must be entered in the user input area, they are 
not calculated by the program. A more detailed description 
of the various outputs is found under the appropriate subrou- 
tine description. The frequency with which output is created 
is controlled by the variable SKIP whi.ch is the number of time 
steps between calls to output routines. 
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2 . 



The "TIME" Routine 



The maximum allowable time step is governed by the 
CFL condition. Simply stated, this means that the time step 
must be small enough so that the characteristic trajectories 
remain within one spatial interval during the time interval. 
The minimum time step is calculated by computing 

Del t = H/ABS [Q+A] 

at every node, where H is the spatial interval, and selecting 
the minimum value of Delt. 

3 . The "TRAK" Routine 

Shock locations at the unknown time level are deter- 
mined by computing the shock speed at the known time level, 
as outlined in Section 2, and multiplying by the time inter- 
val computed by "TIME." The time step is reduced, if neces- 
sary, to limit shock travel to one spatial interval. The node 
immediately to the right of the shock (upstream) is flagged 
for later use by "SWEEP" and "JUMP." 

In calculating the shock speed, the assumption is made 
that the conditions immediately adjacent to the shock are the 
same as the conditions at the nodes to the left and right of 
the shock. 

4 . The "SWEEP" Routine 

The "SWEEP" routine makes the necessary calculations 
to solve Eq. (14) . This involves interpolation at the known 
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time level for the values of Aw, calculation of z, and 
integration of z for each node. 

To calculate Aw, the assumption is made that the 
characteristic trajectories are straight lines. An initial 
guess of the characteristic slope (A) is made and by forcing 
the characteristic to pass through the node at point B in 
Fig. 1, As = A xdelt. The values of q and A at point A are 
found by interpolation and the slope of the characteristic 
is then calculated using the values of q and A at point A, 
This slope is compared with the initial guess and if the two 
are not in agreement, the procedure is iterated until they 
converge. Once As is accurately determined, the values of Aw 
can be calculated by interpolation. Linear interpolation is 
used in the present version of EULER-1. This procedure is 
carried out for each characteristic at each node. 




Figure 6. "RSHOCK" Case 

A deviation from the above procedure is necessary in 
the case where a shock exists in the spatial interval to the 
left or right of the current node. When a shock exists to 
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the right, and the flow is subsonic, as in Fig. 6, care must 
be taken not to interpolate for point 3 based on the slope 
between points 1+1 and I, which would not be accurate due 
to the discontinuity. In such a case "RSHOCK" is called and 
the interpolation is based on the slope between point T and 
1-1. Similarly, when the shock is to the left of the current 
node, "LSHOCK" is called and the interpolation is based on 
the slope between 1+1 and I. 

The calculation of z includes the calculation of 
spatial derivatives of q and A for characteristics 1 and 3 
of Fig. 6. In general, the derivatives associated with 
characteristic 3 are forward differenced and those associated 
with characteristic 1 are backward differenced in keeping 
with the principle of domain of dependence. If the flow is 
supersonic or if a shock exists to the right of the current 
node, all derivatives are backward differenced. If a shock 
exists to the left of the current node, all derivatives are 
forward differenced. Once the derivatives are known, z is 
calculated from Eq . (11) using the average value of A between 

points 1 and I or 3 and I as appropriate. Since the deriva- 
tives are linear, this results in an average value of z 
over the same interval. 

The integration of z is transformed from a time to a 
spatial integration as described in Section II and the 
trapezoidal rule is used to carry out the integration. In 
EULER-1, this has been done in one step using the average 
value of z described above. 
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Equation (14) is then solved by addition of Aw and 
the results of the integration, and the resulting 6w is 
stored until all nodes are calculated. After all nodes have 
been calculated, the variable arrays (w) are updated for the 
next time interval. 

5. The "LSHOCK" Routine 

As discussed above, the "LSHOCK" routine modifies 
the basic EULER-1 procedure at an interior node when a shock 
exists to the left of the node. Interpolation of quantities 
to the left of the node are computed based on the derivatives 
of the quantities to the right of the node. Although this 
assumes that the derivatives do not change between adjacent 
spatial intervals, it is necessary to avoid taking derivatives 
across discontinuities in the flow. 

6. The "RSHOCK" Routine 

Similar to "LSHOCK," "RSHOCK" bases interpolations 
of quantities to the right of the node on the derivatives of 
the quantities to the left of the node when a shock exists 
to the right of the node. 

7 . The "JUMP" Routine 

The "JUMP" routine is used to calculate the conditions 
downstream of a stationary shock as described in Section II. 

If a shock is known to have crossed a spatial node during a 
time interval, which is known once "TRAK" has been called for 
that time interval, the entire "SWEEP" sequence is skipped for 
the node which was crossed by the shock and the conditions at 
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that node are determined using the normal shock relations. Note 
that the node in question is the node downstream of the sta- 



tionary shock. As in "TRACK," it is assumed that the condi- 
tions at the nodes upstream and downstream of the shock are 
the same as the conditions immediately adjacent to the shock. 

8 . The Output Routines 

There are four subroutines which produce various types 
of output as the user desires. "BORDER" and "PLOT" produce 
a graphical presentation of the pressure, density, velocity 
and entropy distributions at selected time levels as seen 
for example, in Fig. 7. "EXACT" produces a plot of the den- 
sity distribution at one selected time interval and compares 
it with the exact solution as shown in Fig. 8. At selected 
time levels, "LIST" produces a tabular listing of the Riemann 
variables, modified entropy, pressure, density and velocity 
distributions, elapsed time, shock speed and discontinuity 
locations. The listing is written to the user's permanent 
storage disk. 

When the value of GRAPHS is set equal to one in the 
"MAIN" routine, "BORDER" is called once to set up the plot 
axis, labels, and headings. "PLOT" is called every SKIP 
time steps to draw the four distribution curves. 

When the value of GRAPHS is set equal to two, "EXACT" 
is called every SKIP time steps to plot the density distri- 
bution compared with the exact solution computed at six 
points of interest. The points are the two end points, which 
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are simply the initial conditions, the head and tail of the 
rarefraction wave, the point just left of the contact surface 
and the point just left of the shock. The spatial locations 
of these points are computed by "EXACT" based on the elapsed 
time and the known values of the wave velocities entered in 
the "MAIN" routine. The exact values of the densities at these 
points must also be entered in the "MAIN" routine. 

When the value of GRAPHS is set equal to zero, the 
tabular listing as described above is sent to the user's 
disk. No graphical output is created, 

D. CAPABILITIES AND LIMITATIONS 

The present version of EULER-] is set up to solve a shock 
tube problem with a single centered diaphragm. Boundary 
conditions for the ends of the tube have not been incorporated 
so the problem must be stopped before the waves reach the end 
of the tube. 

The left side of the diaphragm is the high pressure side. 
The program wil.l not run with the right side as the high 
pressure side without changes to some of the shock correction 
logic . 

The user may select any odd number of grid points limited 
only by the amount of memory available. 

The program tracks shock waves and makes shock jump calcu- 
lations at the appropriate locations but the program can also 
run without the shock tracking feature and jump calculations 
if so desired, with some smearing of the shock discontinuity 
and complete loss of entropy change across the shock. 
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TABLE 1 





LIST OF EDITABLE PARAMETERS 


Parameter 


Function 


N 


Number of grid points 


M 


Controls which discontinuities are 
tracked 

1 = Shocks only 

2 = Shocks and contact surfaces 


GRAPHS 


Controls the form of the output 

0 = Tabular listing 

1 = Pressure, density, velocity, 

entropy plot 

2 = Exact solution comparison for 

density 


SKIP 


Number of time steps between output 
calls 


JSTOP 


Number of time steps calculated 


TRI 


Initial temperature ratio 


PRI 


Initial pressure ratio 


DRI 


Initial density ratio 


QLI 


Initial velocity left of the diaphragm 


QRI 


Initial velocity right of the diaphragm 


G 


Ratio of specific heats 
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IV. 



TEST RESULTS 



EULER-1 was tested on the Riemann shock tube problem. 

The initial conditions and the courseness of the computa- 
tional mesh were varied. Additionally, one run was made 
without the shock tracking and correction features. 

Figures 7 and 8 show the results for initial pressure 
and density ratios of 5, uniform temperature and with the air 
initially at rest, using a grid of 101 points. The exact, 
analytically predicted conditions are also shown in Fig. 8. 

It can be seen that all wave velocities are correctly computed 
and the shock is defined within one interval . The rarefrac- 
tion waves are slightly smeared and the contact surface is 
greatly smeared. A slight transient instability to the right 
of the diaphragm is noticeable, particularly in the plots 
of pressure and velocity. 

Figures 9 and 10 show results for initial pressure and 
density ratios of 1.3. The shock is still well defined in 
the correct interval and the contact surface is not as badly 
smeared as in Fig. 8. There is a loss of accuracy however, 
with respect to the tail of the rarefraction wave. 

Figure 11 shows the results of the original problem with 
a grid of 51 points. As might be expected with a coarser 
mesh, the discontinuities are less sharply defined although 
it is clear that the wave velocities are accurately computed. 
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Figure 12 shows the results for the original problem with 
a grid of 901 points. It can be seen that the EULER-1 
solution closely approximates the exact solution. 

Figures 13 and 14 show the result of running the original 
problem without the shock tracking and correcting features. 
Note that the shock position is incorrectly computed and it 
is smeared slightly. Also note in Fig. 13, the lack of any 
change in modified entropy across the shock. 

Run time for the 101 point mesh is 0.00049 seconds per 
node per time step. The results for Fig. 7 took 50 time 
steps for a total time of 2.46 seconds. For the 901 point 
mesh, the run time decreased to 0.00023 seconds per node per 
time step due to the less frequent use of the shock tracking 
and correction routines per node calculated. The results for 
Fig. 12 took 470 time steps for a total time of 71.85 seconds. 
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SHOCK TUBE RESULTS 

FIRST ORDER N - 101 
DENSITY RRTIO - 5 TEMP RRTIO - 1 
PRESSURE RRTIO = 5 
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Figure 7. 



EULER-1 Results for Pressure Ratio 5 with 
Shock Tracking 
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Figure 8. Exact Solution Comparison for Pressure Ratio 5, with Shock Tracking 
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Figure 9. EULEF-1 Results for Pressure Ratio 1-3 with 
Shock Tracking 
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DENSITY DISTRIBUTION 
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Figure 10. Exact Solution Comparison for Pressure Ratio 
with Shock Tracking 
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Figure 11. Course Mesh Solution for Pressure Ratio 
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Figure 12. Fine Mesh Solution for Pressure Ratio 
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Figure 13. EULER-1 Results for Pressure Ratio 5 
without Shock Tracking 
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without Shock Tracking 



V. DISCUSSION 



A. RESULTS 

The transient instability noted in the results was not 
investigated since the amplitude was small and the effect 
would be of no consequence in most applications. 

At low pressure ratios, the velocity of propagation of 
the tail of the rarefraction wave approached sound speed 
and, in these cases, the EULER- 1 code did not accurately 
compute its location. The reason for this was not investi- 
gated here; however, due to the low pressure ratios at which 
wave rotors can operate > this is of interest and should be 
investigated in future work. 



B. SPECIAL CONSIDERATIONS 
1 . Modified Entropy 

In the QAZlD method, and subsequently, in the EULER-1 
code, the variable S, which here has been referred to as the 
"modified entropy" and which is referred to in Reference 6 as 
simply the "entropy," does not behave thermodynamically as 
one would expect entropy to behave. As a fluid crosses a 
shock wave, the modified entropy of the fluid decreases. 

This is the result of the definition of Eq. (3), repeated here, 
that 



dS 



1 

Y 



dQ R 

T 
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for this varia- 



The use of the word "entropy" and symbol "S" 
ble is confusing since the sign of the change in "S" is in 
direct conflict with well-established conventions in thermo- 
dynamics. For example, Corollary 6 to the 2 n< ^ Law of 
Thermodynamics [Ref. 8:p. 88] would have to be changed to 
read, "The 'modified entropy' of an isolated system decreases 
or in the limit remains constant." At a minimum perhaps, the 
symbol 5 to denote the "modified (negative and scaled) entropy" 
would be helpful. In practice of course, the change in the 
conventionally defined entropy can be recovered from the change 
in the modified entropy, simply by multiplying by (-y) . 

2 . Moretti's Methods 

Incorporating Moretti's methods of handling discon- 
tinuities [Ref. 7] into the EULER-1 code required special 
care . 

There is a fundamental difference in procedure in 
that the EULER-1 code is based on the high pressure side 
being on the left whereas Moretti's formulations are all 
based on the high pressure side being on the right. Although 
essentially a matter of bookkeeping, it is an area of poten- 
tial confusion. 

Another difference is that the Riemann variables used 
by Morreti. are not the Riemann variables used in the QAZ1D 
method. This can lead to difficulty. In fact, in the 
application of Moretti's shock tracking scheme to the EULER-1 
code, the pressures and densities have to be non-dimensional ized 
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by their values on the low pressure side of the shock in order 
to obtain Eq. (14) , and for the procedure to work. 

C. SUITABILITY FOR WAVE ROTOR APPLICATIONS 

Further work is necessary before the suitability of the 
QAZ1D method for wave rotor applications can be fully evalu- 
ated. The following areas need to be addressed. 

1 . Boundary Conditions 

In the EULER-1 code, boundary conditions at the ends 
of the passage have not been incorporated. The code calculates 
the changes of the interior nodes only, skipping the two end 
nodes. In particular, solid wall boundary conditions and 
open-end boundary conditions must be incorporated. No particu- 
lar problems are expected in the boundary conditions themselves. 
However, additional logic may be necessary in the handling of 
discontinuities, such as the reversal of conditions to the 
left and right of a discontinuity after reflection from a 
boundary. The additional logic, which may be considerable, 
is warranted by the simplicity of the QAZ1D method. 

2 . Contact Discontinuity 

No special attention is given to the contact discon- 
tinuity in the EULER-1 code and hence the discontinuity is 
smeared. In a wave rotor application, this discontinuity 
will have to be sharp. Moretti's methods again appear to be 
applicable here and the additional logic needs only to be 
formulated and incorporated into the EULER-1 code. 
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3 . Quasi-One-Dimensional Modeling 



EULER-1 is a one-dimensional code. The solution of 
flows in passages of varying area may be desirable for some 
wave rotor applications. Such problems can be solved in a 
quasi-one-dimensional manner by the addition of the appro- 
priate area variation term to the equations , and the need to 
solve the fully two-dimensional equations is avoided. The 
area variation must be incorporated into the EULER-1 code and 
the code tested on quasi-one-dimensional problems. 

4 . Other Potential Extensions 

Another capability which can be incorporated is the 
ability to handle two gases with different specific heat 
ratios. Eventually, the code should be extended also to a 
second order, one-dimensional version and then to a first 
order, two-dimensional version. 
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VI. CONCLUSIONS 



The development of the EULER equations in the natural 
streamline coordinates using extended Riemann variables was 
reviewed in detail. The QAZ1D solution method, with the 
addition of shock tracking features, was implemented in a 
first order, one-dimensional FORTRAN code with the intent of 
evaluating the method's suitability for wave rotor applica- 
tions. The code (EULER-1) was tested on the shock tube 
problem, with good results. The incorporat j on of boundary 
conditions, an improvement in the contact discontinuity 
definition and the addition of an area variation term for 
quasi-one-dimensional modeling are considered to be neces- 
sary before the QAZ1D method can be accepted as being suitable 
for wave rotor applications. The additional logic required 
for these extensions may be considerable but the development 
is warranted by the overall simplicity of the QAZ1D method, 
and its straightforward extension to viscous, multi-dimensional 
flow modeling. 
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APPENDIX A 



DERIVATION OF THE GOVERNING EQUATIONS FOR 
THREE-DIMENSIONAL INVISCID FLOW (WITH AREA CHANGE) 



The governing equations of a 3-D, inviscid, compressible, 
unsteady flow are derived here in a natural streamline 
coordinate system. The equations are then recast along 
characteristic trajectories in the (s,t) plane and expressed 
in extended Riemann variables, reducing the system of par- 
tial differential equations to a system of ordinary differ- 
ential equations. 



A. DEFINITION OF VARIABLES 
A 



A 

C 

P 

c 

v 

da 



s , n,m 



R 



q 

r 



Cross-sectional area of a differential stream 
tube 

Speed of sound 

Specific heat at constant pressure 

Specific heat at constant volume 

Normal vector for differential area of control 
surface 

Specific internal energy 

Unit vectors in the s f n, and m directions 

Static pressure 

Modified Riemann variable 

Reversible heat transferred 

Velocity magnitude 

Position vector 
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R Modified Riemann variable 

R,, Gas constant 

S Modified entropy 

T Static temperature 

t Time 

u Velocity relative to a standing shock wave 

v Specific volume 

vol Differential element of volume 

V Velocity vector 

W Incoming Mach number relative to a standing 

shock wave 

p Density 

y Ratio of specific heats 

d , (p Flow angles with respect to reference coordinate 

planes . 

V Vector operator 
B, THE COORDINATE SYSTEM 

The natural streamline coordinate system (s,n,m) is 
shown with respect to a fixed rectangular cartesian system 
(x,y,z) in Fig. Al . The system is a curvilinearlv translating, 
right hand orthogonal system which translates with a fluid 
particle along a streamline such that the 's' coordinate is 
measured in the direction of the flow. The ' n ' coordinate 
direction always lies in the plane defined by the 's' coor- 
dinate direction and the fixed ' y' axis. The 'm' direction is 
norma] to the (s,n) coordinate surface. The flow angles, 

0 and (p , are defined as shown in Fig. Al . 
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Figure Al . The Natural Coordinate System 

C. VECTOR OPERATORS IN THE NATURAL COORDINATE SYSTEM 

1 . ?( ) 

If dr is the position vector of a fluid particle 
dr • V ( ) = d ( ) 

In natural coordinates, Eq. (Al) becomes 



fids + idn + i dm]*[( )i + ( )i + ( )i ] 
•s n m 1 s n m 



S! Us + illdn + A<-Um 



3 s 



3n 



?m 



By inspection, 



(Al) 
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V( ) 



(A2) 



a ( ) v , ILli + 1L_L i 

3s s 3n n 3m m 



2. W( ) 



Since 



using Eq . (A2) 



V = q i 



V- V ( ) = 



3 ( 



3s 



(A3) 



3. V-V 



From Eq. (A2) with V = qi 



• V7 = M 



V-V = 



3s 



( A4 



(V-V) V 



-r 



From Eq. (A3) with V = qi 



^ y 

(V-V)V = q £ = 



3(qis ! 

9s 



~ - 9i 

r 9q - , 5 

= C|1 3i 1 s + q 3F 



(as: 



The change in the unit vector i with respect to s is 
derived below and is illustrated in Figs. A2 through A6 . (A 
three-dimensional model is very helpful in visualizing the 
vector geometry.) 

Figure A2 illustrates the natural coordinate system 
translating along a streamline from point A to point B in 
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i +di 
n n 




Figure A2 . Unit Vectors Moving Along A Streamline 




Figure A3. Superposition of i Vectors from Points 
A & B 
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Figure A5 . X-Z Plane 



Y 




Figure A4 . Plane Formed by B Vector and Y Axis 




Plane Formed by a Vector and Y Axis 



h 2 



F ' gure Ab . 



the inertial cartesian coordinate system. Figure A3 is the 

A 

superposition of the i unit vector from points A and B of 
Fig. A2 . The intent here is to express the change in the 
unit vector i s in terms of components along the original 
(s,n,m) directions shown at point A in Fig. A3. Vector OA 

A 

in Fig. A3 is the original i direction and vector OB is 
3i 

3 

i + — — ds . Points C and E are the projection of points B 
s 3 s 

and A in the (x,z) plane respectively. Point D is the projec- 
tion of point C onto the line OE . 

Figure A4 is the plane formed by the OB vector and the 
Y axis. The length of OB is unity by definition. The length 
of BC is sin(0+d6) . 

sin(6+d9) = sin 0 cos dfl + cos 0 sin d0 

which, since dB is a small angle, is 

sin(0+d0) = sin R + cos 0 d0 (A6) 

The length of OC is cos(fl+d0). 

cos(8+d0) = cos 0 cos dfi - sin 0 sin d0 

which, since d0 is a small angle, is 

cos(0+dfi) = cos 0 - sin 0 d6 ( A 7 ) 
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Figure A5 is the (x,z) plane in which the angle cj> is 

/\ 

defined. The i direction is shown at point C. The length 
of CD is OC sindi> which, with Eq. (A7) is 

[cos 0 - sin 0 d0]sin dtp 

which, since dcf> is a small angle, is 

cos 0 dtj) - sin 0 d0 d<£ 

Dropping the second order term 

CD = cos 0 dtp (A8) 

The length of OD is OC cos di[) which, since dtp is a small angle, 
is OC. 

Figure A6 is the plane formed by the OA vector and 

the y axis. GD is the projection of BC and since both are 

vertical lines, BC and GD have the same length, which has 

been determined in Eq. (A6) . The i and i directions are 

s n 

shown at point A. The lengths to be determined are AF and GF , 
since these are the components of AB in the s and n directions. 
For small angles, OD was shown to be equal to OC and GD = BC, 
so that the angle GOD is equal to angle BOC = 0+d0. Angle GOF 
is therefore d0, the length of OG is unity and the length 
of GF is sin d.0, which for small angles is d0 , and it is in 

A 

the i direction. The length of AF is OA -OF. OA is unity 
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by definition and OF is cos d0, therefore AF, the component of 
AB in the i direction is equal to 1 -cos d0, which, since 
d9 is a small angle, is zero. 

The component of AB in the i direction is GF or 

/N 

d0 i n (A9) 

/s 

The remaining component of AB in the i^ direction is 

CD or 



cos 0 dtp i 



m 



(A10) 



Therefore, the vector AB can be written as 



3i 

c 

9s 



ds 



d0 i + cos 0 dd> i 
m T m 



90 , ? 

= — ds 1 + cos 0 -rf- ds 1 

9s n 9s m 



and, dividing by ds, finally 



9i 

ai 



9 0 i Q 3* 7 

7T~ i + cos 0 1 

9s n 9s m 



(All 



Note that Fig. A3 can be considered also to describe 

a change in the unit vector i at a fixed location, with 

respect to time. In this case, the vector AB is equal to 

9 1 

dt and with Eq . (A9) and Eq. (A10) 
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(A12) 



3 i 

c 



3 0 2 Q 3<f) 2 

-r-r l + cos 0 i 

o t n 3 1 m 



With Eq. (All), Eq. (A5) becomes 



(VV)V = 



3q 2 ^ 2 30 2 2 0 3* 2 

q i_ + q l + q cos 0 ■r- L i 

^ 3s s ^ 3s n ^ 3s m 



(A13) 



D. CONSERVATION OF MASS 

In the quasi-one-dimensional differential stream tube, 
shown in Figure A7 , where i is always in the direction of 
V, p, q, and the cross sectional area A are known at the center 
of the element and p and q are assumed constant on any given 
cross section. 



P - 



3p ds 
3s 2 



q - 

A - 



3q ds 
3s ~2 
3A ds 
3s 2 




r\ 


x 


3p 


ds 


P 


i 


3s 


2 






3q 


ds 


q 


+ 


3s 


2 


A 


x 


3A 


ds 




T 


3s 


2 



Figure A7 . 



Differential Stream Tube of Variable Cross 
Section 



The statement of cont.inui.ty is 



The change in the mass 
within the control volume 
with respect to time 



The net influx of 
mass across the 
control surface 
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Since the volume is constant, the left side becomes 



A ds 



l£ 

3t 



For the stream tube, mass enters and leaves only from the 
ends so that the right hand side can be written as 



[p - 



3p ds, r 3q ds, ra 3A ds 

3s 2 ' q 3s 2 J 1 3s 2 



r , 3p ds, , , 3g ds, r ^ , 3 A ds- 

tp+ 3?— 1Iq+ 3#-2" !1A+ ^— ■ 



which on expansion becomes 



- pq 



3 A 
3s 



ds - qA 




3q , 

- p A at ds 



3p 3q 3A , 3 
3s 3s 3s 



On dropping the higher order terms and combining with 
the left hand side. 



iR 

3t 



+ 



pq 



I Ih- 

A 3s 



+ 



3p 

q 3i 




0 



( Ai 4 ) 



Equation (A14) is the form of the continuity equation 
which is required to model a flow as being one-dimensional 
with area change. In general, however, ^ must be expressed 
in terms of q, 0 and cf> . 

With reference to Fig. A8 , the change in the cross- 
sectional area of the stream tube can be written as 



3 A 
3s 



ds 



[dm + 



3 dm 
3s 



ds] [dn 



3dn 
3 s 



ds] 



A 



( A15 ) 
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1 




n 





dn < 



A 



l 



m 



dm 

Figure A8 . Cross Section of a Diverging Differential 
Stream Tube 

After expansion, dividing by A and dropping the higher order 
terms, Eg. (A15) becomes 



1 9A 
A 9s 



1 9dn _1_ 9dm 

dn 9s dm 9s 



(A16) 



From Fig. A9 





which, for small angles, becomes 



1 9dn 



9_e 

9n 



( A17 ) 



dn 9 s 



From Fig. A10 





(A18) 
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Figure A9 . (s,n) Plane of Diverging Differential Stream Tube 




Figure A10. (s,m) Plane of Diverging Differential Stream 

Tube 




Figure All. (s,m) Plane and its Projection onto the 
(x,z) Plane 
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With reference to Fig. All, AB = CD, and 



CD 



OD 



OC 



OD sin dm] 
OC 

cos dm] 
ds cos 0 



(A19 ) 



(A20) 



(A21) 



Combining Eqs . (A18) through (A21) and knowing that dm 

is a small angle. 



1 3 dm 

dm 3 s 



cos 




(A22) 



Equation (A16) , with Eqs. (A17) and (A22), becomes 



1 3A 
A 3s 




cos 




(A23) 



With Eq. (A23) , the general form of the continuity equation 
in natural coordinates is 



i£ 4. „ Ip + . 9q 

3 1 



+ q 



3_p 

3 s 



3s 



r 36 

+ + 



COS 0 



111 

3m J 



0 



(A24) 



E. CONSERVATION OF MOMENTUM 

The statement of conservation for an arbitrary control 
volume which is fixed in the reference .frame can be written 

as 



70 



The change in the 




The net influx 




The net force 


momentum of the fluid 




of momentum 


_L 


on the con- 


in the control volume 




across the con- 


1 


trol surface 


with respect to time 




trol surface 







The net 
body 
force 

on the fluid 



Assuming the fluid is inviscid, the only forces acting on the 
control surface are pressure forces. In the absence of any 
body forces (gravity, electromagnetic, etc.), the statement 
becomes 



/// dvol = - // VpV* da - // P da = 0 



(A25) 



With Gauss' Theorem, Eq. (A25) becomes 



/// [ pV] dvol + /// V • [ VpV] dvol + /// VP dvol = 0 



or 



///{-r; -[pV] + V- [VPV] + VP }dvol = 0 



(A26) 



Since Eq. (A26) is true for any volume, the integrand must 
be zero, thus 




V- [vpvi + VP 



0 



(A27) 
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Expanding the first term of Eq. (A27) and using the vector 
identity 

V- [VPV] = V[VpV] + p(V'V)V 



Eq. (17) becomes 



+ v{^ + 
at at 



V-pV} + p(V*V)V + vp = o 



(A28) 



The second term of Eq. (A28) is zero by the conservation of 
mass so that Eq. (A28) reduces to 



p Tj^r + p(V*V)V + VP = 0 



( A 2 9 ) 



Using Eqs . (A2), (A5), (All) and (A12), collecting terms and 

equating each vector component to zero, Eq. (A29) becomes 



Vp 




3q , 3P. 

pq at + V 



0 



\n [ pq 



39 

at 



+ pq 



2 3 6 

3s 



+ 




0 



A 

i 



r 

m ' 



H 



P q COS 0 + 



pq cos 6 



a_£ _ap 

3 s 3m 



0 



(A30) 

(A31) 



(A3 2 ) 



Using the perfect gas relation 

A 2 = I* 

P 
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and the identity 



3P 

3s 



Eqs . (A30) through 



= P 



(A32) 



3 UnP 
3s 

become 



3q 

3 1 



+ 



q 



dd 
3 1 



+ 



q 



3_£ 

3 1 



+ q 



which are in the 



= 

9 s 



36 

3s 



9 cj) 

9 s 



desired 



A^_ 3 inP 

y 3s 

3 tnP 
yq 3n 

A 2 3 jjnP 
yq cos 0 3m 



form. 



(A33) 



(A34) 



(A3 5) 



F. CONSERVATION OF ENERGY 

The conservation of energy for an arbitrary control volume 
is 



The rate of increase of 
energy within the con- 
trol volume with respect 
to time 



The rate at which energy 
is entering the control 
volume across the 
boundary 



The net rate of 
work done on the 
fluid at the 
boundary 



(A3 6) 



Neglecting gravitational potential, the energy per unit mass 
is 
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where e is the internal energy per unit mass and q is the 
velocity magnitude. For an inviscid fluid in the absence of 
body forces, the only work done on the fluid is pressure work. 
Thus Eq. (A36) is written mathematically as 



/// [e P dvol = — / / [ e + ^-J pV-dit - // PV-da 



(A37) 



Using Gauss' Theorem, Eq. (A37) is written as 

2 2 
/ [e + 2— ] p dvol = - fffv • [ (e + ^pV] dvol 

- /// V. [Pv] dvol 

or 

2 2 

-|^{[e+3j-]p} + V-{[e+^-]pv} = -V- [Pv! 

Expanding , 

2 „ 2 2 2 
|^[e + + p |^-[e + —] + ( pv) V • [e + 5— ] + [e + ^— ] [V • ( pV) ] 

= -V- [PV] 
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or 



(A3 8 ) 



[e + \ ] [ ft + V ' (pV) ] 



+ p h [e + V> 



+ (pv) V- [e +^-] 



= -V- [PV1 



The first term of Eq. (A38) is zero according to the conserva- 
tion of mass. Expanding the remaining terms of Eq. (A38) , 



P 




+ 





-V • [PV] 



or 





£ V- [PV] 



(A39) 



Equation (A39) can be written as 



i?® + 

Dt Df 2 J 



-(V* VP) - -(V* V) 

P P 



(A40) 



With V = qi s in the natural coordinate system. Eq . (A29) 

becomes 



3 V 

3 1 



+ 



q 



3V 

3s 



1 

P 



VP 



or 



DV 

Dt 




(A41) 
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Taking the dot product of V with Eq. (A41) , 



- DV 
V " Dt = 



- |[V-VP] 



(A42) 



and substituting Eq. (A42) into Eq. (A40) 



De + D_ f q_ 
Dt Dt ' 2 



77 Dv P . „ — . 
V * Dt ' p (V ‘ V) 



(A43) 



But 



DV 

Dt 



D(qig) 

Dt 



Dq : + Di s 

Dt X S q Dt 



ga i + 

Dt S 



3i 



3i 



qt 3t^ + q 



3 s 



and using Eq , (Ail) and Eq. (A12) , 



DV Dq 7 .30 ? , „ 3<h 7 , 

Dt = 1 B + q f^T + cos 6 at 



Dt s 



3t n 



3t m' 



+ q 2 (|i i + cos 0 M i ] 
^ 3s n 3s m 



Therefore 



- DV 
V *Dt 



f Dq 



qi i 

^ s Dt s 



i + 



, 3 0 2 3 0 '.' 

n +< ! 3 V 1 V. 



2 ff U n + q 003 6 [ lt + q H 1 



= a 



Dq = 
Dt 



s_ t a!] 

Dt 1 2 1 



Therefore Eq. (A43) reduces to 
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De 

Dt 



(A44) 



= - | (V ' V > 



From the conservation of mass 



|| + V- (pV) = 0 



or 



5 p 

St 



+ 



V- Vp 



+ 



P(V-V) 



0 



which, in the natural coordinate system gives. 



l£ 

St 



+ 



q 



i£ 

Ss 



-p(v- V) 



or 



I ££ 

P Dt 



(V - V) 



(A45) 



Substituting Eq . (10) into Eq . (9), 



De P Dp 

Dt 2 Dt 

P 



(A46) 



The definition of the modified entropy is 



dS 



1 

Y 



dQ„ 

R 

~¥~ 



11 



and the first law of thermodynamics for an elemental mass 
requires that 

dQ = de + Pd (— ) 

R p 

so that 



dS 



— ^ [de + Pd (-) ] 
yT p ' 



= i[de + 

yT 



-7 dp] 

P 



and Eq. (A46) implies that 



DS 

Dt 



0 



(A47) 



which states that modified entropy is conserved along streamlines 

G. MODIFIED ENTROPY EVALUATION 

With the modified entropy defined as 



or 



dS 



YT 



S B " S 



A 



B dQ 

/ 

A 



(A48) 



From the first law of thermodynamics 
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de + P dv 



(A49) 



d Q R = 



where e is the internal energy and v is the specific 
Substituting Eq. (A49) into Eq . (A48), 

„ o _ 1 f B de + P dv _ 1 r r B de A f B P dv, 

a " s b - v a j T ■ y' a j t h > 

For a perfect gas 

P v - R„T and de = C dT 

Cj V 



for which Eq. (A50) becomes 



q _ q 
A B 



B C v ^dT 

Y Y T 



B R_ dv 

/ A-] 



v 



which after integration yields 



*1 T v 

S* = S D + -[C £n =^- + R In —} 

A B y L v T a v a ' 



Substituting 



R 



G 



C 

v 



(Y-l) 



into Eq , (A51) 



volume . 



(A50) 



(A51) 
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T v v 

s + £[C Un -1-^n -®) + C v y£n -A 

B Y v T a v ft v a 



S + -[— y £n(A — — -=T— — -) + -A £n— ] 

B Y Y-l P A V A R B V B ^ P B 



and since R = R 

g a g b 



S A S B + y(Y-l) [?n P A 



P B 

Y £n— ] 

'A 



or 



R. 



-I + i_ 

R g Y(Y-l) 



p p„ 

[In - in -A.] 



p b y 



P A Y ' 



or 



R 



G 



r g y(y-D 



P P 

[ £n — - - £ n — ] 



p A Y 



p b y ' 



(A52) 



If the reference condition is chosen to be 




2 

Y-l 



1 

Y(Y-l) 




Eq. (A52) becomes 



S 

R 



B 

G 



2 

Y-l 



1 

Y (Y-l) 



in 



B 



P B Y 



(A53) 
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H, FURTHER TRANSFORMATIONS 



Equations (A24) , (A33) through (A35) and (A47) are in the 

required form. However, some of the variables need to be 
transformed, namely, derivatives of An P need to be expressed 
in terms of A and S and derivatives of p must be expressed in 
terms of q and A. 

Using the first law of thermodynamics expressed as 

dQ = dh - — dP 

P 



Eq. (A48) may be written as 



-y T VS 



Vh 




VP 



A 

or for the i component 



3S _ 1 9C p T P_ 3 An P 

3s T 3s pT 3s 



(A54) 



Assuming to be constant and using the equation of state 
for a perfect gas Eq. (A54) becomes 



3S _ ^p 3 r P i _ 3 An P 

3s t 3s R p 3s 



and since C = y r /(Y -Y ) 

P 

3S _ yR G pR G 3 f P i 3 £n P 

Y 3s (y-1) P 3s L R g p J R G 3s 
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Since R is a constant it may be moved in or out of the 
derivatives giving 



- Y 


93 V 


_J_ P iL_ r£i _ 

Y-l P 3s 1 p J 


3 Jin P 


3s 


and since y is also 


constant, similarly 


- Y 


— [— ] = 
3 s R r 


Y o 3 rPy, 

v-1 Py 3s 1 p J 


3 Hn 
3 s 


,2 

Using A 


= YP/P 






“ Y 


— [_§_] 
3s V 


Y 1 3 A 2 


3 Jin P 


y-l a 2 3s 


3 s 






2 Ay 3A 


3 Jin P 






, n N , 2 3s 

(y-i) A 


3s 


and i on 


rearranging 








3 9, n P 

3s 


2y 1 3A 

i -f V 

Y-l A 3s 1 


— [5.1 

3s 1 R' 



Substituting Eq. (A55) into Eq. (A33) yields 



M +a ia + -2A.3A + A 2 i_.s, = 

3t q 3s y- 1 3s 3s l R J 



and the derivative of J,n P has been removed. 



(A55) 



(A56) 
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Since 



A 



2 



yP 

P 



logarithmic differentiation yields 




(A57) 



For a pure substance 



P = P(p,S) 



and an isentropic process 



P = P ( p ) s 



Differentiating 



= 1L - a 2 

9 p S dp 



therefore 



dP 

P 




Substituting Eq. (A58) into Eq. (A57) gives 



(A58) 



2 



dA 

A 



dp_ _ dp_ _ _ d £ 

P P Y “ P 



or 
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dp 



2 p dA 
(Y-l) A 



Along a streamline 



dp 




q 



Ip 

3s 



and since the fluid is the same from streamline to 
Eqs . (A59) and Eq , (A60) can be combined giving 



3_p 

at 



+ 




2 p dA _ 2 p , 9 A 3jA 

y-l ~A ~ (y-l)A at q 3s 



Substituting Eq. (A61) into Eq. (A24) yields 



3A 

at 



3A 

3s 



y-l 



A la 

A 3s 



+ q A 



y-l r 3 6 
2 1 3n 



+ cos 0 



Mi 

3m J 



0 



and the derivatives of p have been removed. 

The governing equations are now summarized: 



iA + a £A + izi A lq + u A izirli + rn , fi 111 
3t q 3s 2 A 3s q A 2 -3n cos 9 3m J 



= 0 



Ml + a Ml + ^a d_A + ^2 a rS, 

3t "3s y-l 3s 3s l R J 



= 0 



36 



30 



at + q as 



A 3 en P 
yu 3n 



11 + a 11 = 
at q ‘ 



0 



3 in P 



yq cos 0 3m 



(A5S ) 
(A60) 

streamline 

(A61 ) 

(A62) 
( A 6 3 ) 

(A64) 

(A65) 

(A66) 
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Note that Eqs . (A65) through (A67) are in the form 



3X 3X 

3 1 3s 



Z 



If X = X (s , t) , then 

... 9X ^ ^ 9X , 
dX = tt dt + ds 

a t dS 



and 



dX _ 3 X ds 3X 

dt 3t dt 3s 



(A68) 



where (ds/dt) is the "character] sti c" direction in the (s,t) 
plane along which Eq. (A68) holds. Along this di rection , the 
equations may be transformed from partial to ordinary differ- 
ential equations. What is desirable is to find transformed 
variables in terms of which Eqs. (A63) and (A64) will be in 
the same characteristic form. For this reason the modified 
Riemann variables 



Q = q + As/Pq 



(A69) 



= q - AS/R g 



R 



(A70) 



are defined, where S/R„ is defined in Eq . (A53) . The modified 

Riemann variables can be introduced by manipulation of the 
equations as follows. First, multiplying Eq. (A63) by S, 
gives 



S + q S fs + AS + q A s + Cos 0 |i] = 



and multiplying Eq. (A67) by A gives 



(A71) 



. 3S . . 3S n 

A 3t + q A 3 ? = 0 



(A72) 



On adding Eqs . (A71) and (A72) to Eq, (A64) , and introducing 



A || - A || = 0 

3 s dS 



(A73) 



and 



* S I? - A * ff - 0 



(A74) 



after appropriate rearrangement 



3 q 

3t 



+ A 



as 

at 



+ s 



3 A 

at 



+ CT 



aq 

as 



, as aA 
+ q A a? + q s ai + A 



aq 

3 s 



+ A‘ 



as 

as 



. „ a a 
+ A s ai 



Y-l 

2 



A 




+ 



, „ 3A 
A s as 



A 



aq 

as 



2 A 3 A 
y-l as 




q A S f t ; — + cos 

an 



e f£l 

3m 



0 
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Rearranging and introducing Eq. (A69) gives 





Y-l 



qA s [ lf +cos 9 H ] 



(A75) 



2 



By subtracting Eqs . (A71) and (A72) from Eq. (A64), again using 

Eq. (A73) and Eq . (A74), and introducing Eq. (A70) , one also 

obtains 



Together with Eqs. (A65) through (A67) , Eqs. (A75) and (A76) 
form the system of equations which describe the isentropic flow 
of an inviscid perfect gas under unsteady, compressible condi- 
tions, and which may be solved as a system of ordinary differ- 
ential equations along characteristic trajectories in the 
space-time domain. 

I. SHOCK JUMP EQUATION 

The analytical expression for the change in the extended 
Riemann variable, Q, across a stationary shock is derived 
below. 

Consider a shock moving with a velocity V into a fluid 

J 1 s 

with velocity q p as depicted in Fig. A12, with the conditions 






(A76) 
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V 



A 

O- 



-O B 



U = q -V 
A s 

A o- 



“b = q B' V s 
-o B 



Figure A12 



Figure A13 



to the left of the shock being denoted by A subscripts and 
conditions to the right by B subscripts. If a velocity of -V 
is imposed on the entire system, the situation becomes one of 
a stationary shock, as depicted in Fig. A13. In both cases, 
the high pressure side is to the left (A side) . Since all 
velocities are defined positive to the right, a positive value 
for the relative mach number is defined as 



W 




'%-V 



B 



(A77) 



The extended Riemann variable, Q, is defined as 



Q = q + AS 



(A78) 



where 
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If all velocities are non-dimensionalized by A R and pressures 
and densities are non-dimensionalized by their values at B, 
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The ratios of pressure, density and sound speed across a 
normal shock are well known functions of Y and the Mach number 
which in this case is W. They are 



^A _ 2y 2 y-l 

P B Y+l W Y+l 



(A81) 
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P B ( Y -Dw 2 +2 
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(A83) 



An expression for (q^ - q ) /A is obtained using mass 

A B o 

conservation . 
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Subtracting one from both sides. 



U A U B 
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U A U B 
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" , * 2 
(Y+l) W 



Substituting Eq. (A77) for u D in Eq . (A84) gives 

.D 



U A- U B 
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_ 2 (1 - w ) 

T~ 

(y+l) w 



A U B _ 2 ( fry 2 - 1 ) 

A fi (y+l) w 



Since = 



q A - v s and U B = q B - v s 



U A U B 



= ^A- V s - ^B + V s 



q A" q B 



Thus Eq . (A85) can be written as 



(A84) 



(A85) 
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Substituting Eqs. (A86) and (A81) through (A83) into Eq . (A80) 

gives 
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APPENDIX B 



USING EULER- 1 ON THE NPS VM/CMS SYSTEM 

A. MEMORY REQUIREMENTS 

Storage should be defined as one mega-byte to ensure ade- 
quate memory for DISSPIA requirements. Also ensure that room 
exists on the user's disk for the listing file if that 
output option is selected. 

B . TERMINAL REQUIREMENTS 

The type of terminal required to run EULER-1 depends on 
the output option selected. When using the tabular listing 
output, any terminal tied into the VM/CMS system may be 
used. Graphical output may also be created using any terminal 
and the graph stored on the user's disk in a metafile for later 
viewing at a graphics terminal or for printing. To have the 
graph stored on the disk, use the "COMPRS" command on line 
223 of the program and comment out the "TEK618" command on 
line 224. When graphical output is selected and it is desired 
to view the graph at the terminal at execution time, an IBM 
3277-Tek618 dual screen terminal must be used. To use this 
option, use the "TEK618" command on line 224 and comment out 
the "COMPRS" command on line 223. When using this option, a 
low quality hard copy of the TEK618 display may be obtained 
using the TEK4631 hard copy unit attached. 
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C. PROCEDURE AND COMMANDS 

After loging on to VM/CMS , use the command "DEFINE 
STORAGE 1M" to increase the virtual memory as recommended 
above. This will remove you from the CMS environment. Re- 
enter CMS with the command "I CMS." 

Edit the user input area of EULER-1 FORTRAN as desired 
for the particular problem being run. When graphical output 
is selected, comment out either the "COMPRS" or the "TEK618" 
command as desired. 

Compute the program using the command "FORTVS EULER-1." 
After compilation, an EULER-1 TEXT and an EULER-1 LISTING 
file will reside on the user's disk. At execution, if a 
tabular output has been selected, it will be sent to the 
user's disk as "FILE FT09F001." To give this file a particu- 
lar name, say "EULER-1 LISTING," use the file definition 
command 

"F I 9 DISK EULER-1 LISTING A ( PERM" 

at compile time. An alternate and convenient means of com- 
piling the program is to set up an exec file on the user's 
disk by the name "EULER-1 EXEC" which contains the commands 

F I 9 DISK EULER-1 LISTING A ( PERM 

FORTVS EULER- 1 

Then to compile the program and define the output file, use 
the command " EULER- 1." 

After compilation, to execute the program, use the 
command "DTSSPLA EULER-1." The message 
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...YOUR FORTRAN PROGRAM IS NOW BEING LOADED... 
...EXECUTION WILL SOON FOLLOW... 

should appear, followed shortly by 

...EXECUTION BEGINS... 

If tabular output was selected, proper termination will 
result in a ready message. If graphical output was selected 
using the TEK618 option, the plot will appear on the TEK618 
terminal, the program will pause, and a 

"...press ENTER to continue" 

message will be displayed on the 3277 terminal. At this 
point the hard copy must be made, if desired, using the 
TEK4631 hard copy unit. After pressing the ENTER key on the 
3277 terminal, the plot will be erased and the program will 
terminate. Proper termination will be indicated by the message 

"END OF DISSPLAY 9.2 #### VECTORS IN 1 PLOT..." 

and a ready message. 



94 



APPENDIX C 



EULER-1 FORTRAN CODE 
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